Strains and pseudo-magnetic fields in circular graphene rings 
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We demonstrate that circular graphene ring under a shear stress displays strong pseudo-magnetic 
fields. We calculate the pseudo-magnetic field both from continuum elasticity theory as well as 
molecular dynamics simulations. Stable wrinkles are induced by shear deformations and lead to 
enhancement of the pseudo-magnetic field. The strong pseudo-magnetic field found here can be 
observed by imaging graphene flake at the atomic level e.g. through scanning tunneling microscope. 
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I. INTRODUCTION 



Graphene has recently attracted intensive interest as a promising candidate material for the new generation of 
electronics and spintronicsi. One of the exciting physics on graphene is strain exerted on graphene samples — —. It 
was proposed that strain can be utilized to generate various basic elements for all-graphene electronics^. 

Semiconductor quantum rings have been investigated by a number of groups^—. In graphene ring, the spectrum 
reveals signatures of effective time-reversal symmetry breaking, in which the spectra are most naturally interpreted 
in terms of effective magnetic flux contained in the ring, even when no real flux is present. Quantum rings can 
be considered as prototypical devices in mesoscopic physics, as they show one of the most basic coherence effects, 
namely the Aharonov-Bohm effect, oscillations of the transmission as a function of the magnetic flux through the 
ring. The reason for these oscillations is the phase difference between electrons traveling along the different arms 

S of the ring. Furthermore, Benjamin and Pachos 10 proposed creating a ring of single-layer graphene in which they 
induce d-wave superconductivity via the proximity effect or directly make graphene superconducting by doping. The 
quantum qubits would be built around the 7r-junction that naturally occurs in graphene and would not require bilayer 
structures. The Aharonov Bohm oscillations, on the other hand, has been observed 9 - in a graphene ring, consisting of 
a planar honeycomb lattice of carbon atoms in a ring shape, by changing the voltage applied to the side gate or the 
i^i , back gate. 

When the graphene sheet is under tension, the side contacts induce a long-range elastic deformation which acts as 
CN ! a pseudo-magnetic field for its massless charge camer a 11 ' 12 . This is because strain changes the bonds length between 
k*" ■ atoms and affects the way electrons move among them. The pseudo-magnetic field would reveal itself through its effects 
on electron orbits. The tension can be generated either by the electrostatic force of the underlying gatei£ by interaction 
of graphene with the side walls^, as a result of thermal expansion^ or by quench height fluctuations^. A particular 
strain geometry in graphene could lead to a uniform pseudo-magnetic field and might open up interesting applications 
in graphene nano-electronics with real magnetic fields^. Mechanical strains can introduce new environments in 
studying novel physics of graphene. 

It is commonly believed that strains have important influence on the electronic structure of graphene 18 . A graphene 
ring is a particularly convenient geometry. Strains can possibly be manipulated efficiently in samples with good 
adhesion to the substrate, such as graphene layers grown epitaxially on SiC. Recently, the physical properties of 
graphene when its hexagonal lattice is stretched out of equilibrium have been investigated by many group a 19 ' 20 . 
Scanning tunneling microscopy studies on graphene surface have indeed revealed a correlation between local strain 
and tunneling conductance 20 . Motivated by experiments pointing to a remarkable stability of graphene with large 
strains, we have carried out theoretical analysis and the molecular dynamics simulation to explore the pseudo-magnetic 
field in strained graphene ring. 

In this papers we focus on a particular aspect of the physics of graphene ring namely the appearance of gauge fields 
and its corresponding pseudo-magnetic fields which arise in the circular graphene ring when a shear force is applied 
to its boundary. The common belief is that the morphologies of graphene surface under strain are negligible. The 
aim of this paper is to show that, contrary to these expectations, the structure of a deformed surface can lead to a 
strong pseudo-magnetic field. 

The paper is organized as follows. In Sec. II we introduce our model and formalism. In Sec. Ill, our numerical 
results for the strain and the pseudo-magnetic field for deformed graphene ring are presented. Finally, we conclude 
in Sec. IV with a brief summary. 
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II. MODEL AND THEORY 



We analyze both analytically and numerically the strains and pseudo-magnetic fields in circular graphene ring under 
a shear stress. A shear force is applied to the boundary, inducing shear deformations inside graphene ring. 

The properties of graphene presented lie on the special character of its low energy excitations obeying a two 
dimensional massless Dirac equation. The graphene ring with valley degree of freedom, r = ±1 for the inequivalent 
K and K' valleys, is modelled by the massless Dirac Hamiltonian in the continuum mode l 21 ' 22 

Ho = hv-pT (a i k\ + cr 2 fc 2 ) 

where fcj is an envelope function momentum operator, vp ~ 10 6 m/s is the Fermi velocity, and <7i = x,y, z are the Pauli 
matrices that act on the sublattice pseudospin degree of freedom. Hence it is important to establish the robustness 
of the low energy description under small lattice deformations. The concepts of gauge fields and covariant derivatives 
can be translated into the language of differential geometry, based on differential forms. When ideal graphene is 
distorted, the effective Hamiltonian will be changed into 

H = v F ( p-e A) • a (1) 

where the induced vector potential field is defined through the deformations of sample and p is the momentum in the 
polar coordinate. The induced gauge fields can be calculated through the following expressions^ 

v F eA x = hg 2 {u xx ~u vy ) (2) 
vpeAy = —2hg2U xy 

with 



.92 = — t (3) 

K = 

2B 

where t ~ 2.7 eV is the nearest-neighbor hopping parameter and (3 — <91n(i)/<91n(a) ~ 2 is the electron Griineisen 
parameter. For the shear \l and bulk moduli, B we have used 2 ^ the values fx = 9.95 eV A' 2 and B = 12.52 eV A' 2 . 
We thus find that k « 0.56. Once the induced gauge filed is obtained, the pseudo-magnetic field can be calculated by 
B = VxA 

When the graphene is deformed due to the force exerts on the boundaries, the strain tensor can be calculated from 

u a p = - + ^- (4) 

where u is the atomic displacement field and h is the out of plane displacement. The contribution of the out-of-plane 
atomic displacements are noticeable in our numerical calculations. 



A. Analytical expression for pseudo-magnetic filed 



We consider the case of a graphene ring of radius R and width W = B — R\ where the clamped outer and 
inner boundaries are circles of radius R and R\, respectively. Additionally, we assume that the out-of-plane atomic 
displacements being zero, h = 0. The displacement at the outer boundary is ug(R,0) — Ug,u r (R,6) = 0, which the 
shear deformation at the boundary is defined by U$, and at the inner boundary we have ug{R\,9) = u r (R\,6) = 0. 
The displacements in the graphene ring are 



u r (r,6) = 
ug(r, 6) = Ue 



r\r 



Rr 



r(i? 2 -i? 2 ) R 2 -R 2 



(5) 



The second term in the expression for ug{r,9) is a pure rotation, required to satisfy the boundary conditions. The 
only non-zero component of the strain tensor is u r g. The strain at radius r is 



u r g(r, 9) 



1 / dug 1 du r Ug\ UgR\R 



2 V dr + r 89 



r 2 {R 2 -R 2 ) 



(6) 
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The maximum strain at the boundary is actually given by u — UgRj (R 2 — R\). Using polar coordinates (r, 9), Eq. ([2]) 
can be rewritten as 



A r = $ 
+ 

+ 

where $o = h/2e is the quantum unit of magnetic flux and c = v/3k/(2tt) is a constant. Furthermore, the induced 
gauge field in the polar coordinate is given by 



cB du r U r I dug 

— (-5 "St) 00 ^ 3 * 9 ) 

a or r r ou 

dug Ug I dug 

( -¥ + 7-;¥ )smW1 ' 

^[(-^ + ^--^)cos(3^) 
a dr r r dO 

du r u r I dug . 

( T + 7n¥ )smffl ' (7) 



cBU 8 2R4R . 



c c/3C/ e 2i? 2 i? 



Eventually, the pseudo-magnetic field acting on the electrons is simply given by 



It would be worthwhile mentioning that the pseudo- magnetic field diverges near the clamped sites as B(r) 
8c$ PuR 2 /(ar 3 ). 



B. Molecular dynamics simulation 



We used Molecular Dynamics simulation (MD) to simulate a suspended circular graphene ring, consists of a planar 
honeycomb lattice of carbon atoms in a ring shape. Rigidly clamped boundary conditions were employed. We 
simulated the system at different temperatures by employing Nose-Hoover thermostat to help the system reach 
equilibrium at a given temperature. Our present results are limited to T = 50K. In this work, we used both the 
Brenner's bond-order potential 24 i 26 incorporating the second nearest neighbors interactions through their dependence 
on the bond angles and the second-generation reactive empirical bond-order (REBO) potential 2 ^ for the carbon- 
carbon interactions (for more details see Appendix A). It is well established that MD simulations predict the correct 
mechanical properties and the form of the structures with carbon atoms by using the Brenner potential. 

The number of carbon atoms in our simulations is 5125 corresponding to graphene ring of radius R = 7 nm and 
width W = 4 nm. At the beginning of calculations, we simulate the circular graphene ring at a given temperature 
with circular clamped boundary conditions. After reaching a stable configuration, we rotate the outer boundary of 
graphene ring about #0 = 12° as the system shown in Fig. [T| left panel). For later purposes, we name the system as 
case one. In this case, the displacement of the boundary atoms is about Ug = Rtt9q/180. For the sake of comparison, 
we would like to consider another system in which the inner boundary is pulled down. In the later case, in order to 
avoid from strong wrinkles, we reduce the rotation of the outer boundary about 6° and then pull slowly down the 
central part of graphene sheet about z = — 1.7nm as shown in Fig. [TJ right panel). Again, the latter system is labeled 
as case two. 



III. NUMERICAL RESULTS 



From the analytical calculations, the pseudo-magnetic field obtained in the graphene ring geometry where h = is 
shown in Fig. [2] For given the circular radius, B = 7.89 x 10 7 cos(3#) / rr, 3 in units of tesla where rr, scaled in units of 
Angstrom. The maximum field occurs in the vicinity of the inner boundary and it decreases like 1/r 3 . Although the 
analytical model described above can explain the behavior of the strains and magnetic field, we point out that the 
wrinkle structures represent an important feature of the results. 
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We use two different systems in the numerical simulations. We first study a graphene ring which is placed in the 
x — y plane and set at a given temperature by allowing the surface fluctuations in the third dimension, and then 
the outer boundary of graphene ring is rotated while the inner boundary is clamped as is shown in Fig. 1 ( left 
panel). In case one, we have seen noticeable wrinkle structures around the inner boundary along the perpendicular 
direction (z-axis). In case two, on the other hand, the graphene ring surface is smoother, but the average carbon- 
carbon distance is longer than the one obtained in case one. In the following, we will describe two cases and find the 
pseudo-magnetic fields emerge due to the strains. In our MD simulations, the equations of motion are integrated by 
the Verlet algorithm with a time step of 0.5 fs. Our numerical results show that the formation of the wrinkles are 
quantitatively sensitive to the details of the potential. We simulate a system with two mentioned potentials and find 
that the order of the pseudo- magnetic fields are the same ( for more details, see Appendix A ). Therefore, we mainly 
use the Brenner potential^ to calculate the pseudo-magnetic filed. 

A. Deformed structure 

As it can be seen from the left panel of Fig [TJ the main part of the deformation occurs in the area near the inner 
boundary and makes noticeable fluctuations along the z-direction. Due to the appearance of wrinkles in the system, 
the strain increases near the inner boundary and reduces in the far regions. In contrast to the analytical calculations, 
the strain does not behave as r~ 2 . The analytical model does not describe well a system with out of plane wrinkles. 

In the right panel of Fig[TJ we show case two, where the static wrinkles are approximately washed out and seemingly 
the surface of graphene ring is much smoother than case one. The structure of surface displacements can also be seen 
in Fig |31 This figure clearly shows that the morphology of the surfaces are enormously different and we expect to 
have a stronger pseudo-magnetic fields for case one with respect to case two. 




FIG. 1: ( Color online) A representative atomic configuration in MD simulations for circular graphene ring at T = 50 K. Left 
panel: case one in which the rotation of the inner boundary is 12° in x — y plane. Right panel, case two in which the inner 
boundary has been rotated about 6° and afterwards the central part has been pulled down about z — —1.7 nm . 

It should be noticed that the structure of the simulated samples weakly depend on temperature up to T — 300K. 

B. Pseudo-magnetic field 

Since the structure of wrinkles, after the relaxation of graphene's shape, is static, we can calculate the gauge fields 
and also the pseudo-magnetic fields from Eq. ((2J. 

Our numerical results for the induced gauge field calculated from the strains are illustrated in Fig. U for two 
mentioned cases in the x — y plane. We take a time average over the atom positions in order to pass over the thermal 
fluctuation effects and study only the wrinkling structures. The length value of the vector denotes the absolute value 
of the gauge field. For case one, the induced gauge fields are not uniform in the area which wrinkles appear. 

The pseudo-magnetic field can be calculated from the induced gauge field and its distributions are demonstrated in 
Fig. [5] It is worth mentioning that the non-zero pseudo-magnetic field, B z almost occurs along the wrinkles, however 
its value and structure depends on each wrinkle structure. This is due to the fact that the structure of the wrinkles 
is different along different lattice directions such as armchair and zigzag. The value of the pseudo-magnetic field that 
increases with increasing curvature of the wrinkles mostly occur around the central part. 

The structure of the pseudo-magnetic field along one of wrinkles, which is indicated in the inset figure, is illustrated 
in Fig [B] (left panel) for case one. Notice that B z decreases for x < along the wrinkle. The curve can be fitted quite 
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well with the expression exp(— ax) where a = 1.2 ± 0.1, which displays fast decaying of the pseudo-magnetic held 
close to the outer boundary. The pseudo- magnetic held behaves randomly for x > 0. The latter behavior is due to the 
fact that the path is no longer along a specific wrinkle. As it can be observed from the results, the pseudo-magnetic 
field value is immense. Similar massive pseudo-magnetic fields raised by highly strained nanobubbles that form when 
graphene is grown on a platinum surface have been measured by Levy et al££. 

For case two, which wrinkles are disappeared, the symmetry of the graphene is clearly induced in the B z , (see Fig[5j 
right panel ). 

In order to have a better understanding regarding the structure of the pseudo-magnetic fields, we calculate the 
pseudo-magnetic fields for two cases along the azimuthal angel at given radius, r. Fig. [Jj shows the pseudo-magnetic 
field B z as a function of 6 around a ring for radius r = 4 nm in case one (left panel) and moreover the pseudo-magnetic 
field is shown at r = 6.5 nm for case two (right panel). In the latter case, we have seen the behavior most likely as 
B z oc cos(36*). In this figure, we have plotted function, which is proportional to cos(36*) as solid lines. Apparently the 
result in the left panel refers to case one behaves similarly to case two with respect to the azimuthal variable, however 
there are detectable fluctuations due to the appearance of the wrinkles. 

The analytical expression of the pseudo-magnetic field acting on the electrons is B(r, 6) oc \ cos (3$) when atomic 
out-of-plane displacements are ignored. The expression was obtained for a system that lies in the x — y plane, h = 0, 
where there are no wrinkle structures on the flake. Since case two, which has less wrinkle structures, is somehow 
similar to the aforementioned system, we have found numerically the cos(3#) dependence for the pseudo-magnetic 
field. However due to the fact that the inner boundary is pulled down, the pseudo-magnetic held behaviors differently 
as a function of r. 



In summary, we have investigated the strains and pseudo-magnetic fields in circular graphene ring under a shear 
stress. We find, from the elastic theory, the induced gauge filed as function of the maximum strain at the boundary, 
u. The magnitude of the pseudo-magnetic field near the boundary is B m 4v / 3 ( I > o/3 K u/(7rai?). Moreover, The field 
diverges near the clamped sites as B(r) oc 1/r 3 , where r is the distance to the site while the strains diverge as 1/r 2 . 
From numerical simulation results, we find wrinkles structures on graphene flakes and the pseudo-magnetic field is the 
same order of magnitude that obtained by analytical calculations. We also find that the wrinkle structures represent 
an important feature of the pseudo-magnetic field. In addition, it is also shown that the pseudo-magnetic field behaves 
like as cos (30). These results are essential for understanding the electronic properties of graphene ring and its strain 
engineering for potential applications. 



We thank A. Fognini and A. Naji for their useful comments. This research was supported in part by the Project 
of Knowledge Innovation Program (PKIP) of Chinese Academy of Sciences, Grant No. KJCX2.YW.W10. F. G. is 
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We used both the Brenner's bond-order potential 24 incorporating the second nearest neighbors interactions through 
their dependence on the bond angles and the second-generation reactive empirical bond-order (REBO) potential 2 ^ 
for the carbon-carbon interactions. The latter potential is based on the empirical bond-order formalism and allows 
for covalent bond binding and breaking with associated changes in the atomic hybridization. Consequently, such a 
classical potential allows to model complex chemistry in large many-atom systems. The Brenner bond-order potential 
can be written in the following general form for the binding energy, 



The first term is repulsive, and the second one is attractive, rij is the distance between pairs of nearest-neighbor 
atoms i and j. Although this expression is a simple sum over bond energies, it is not a pair potential since the b^, 
which is called the bond-order factor, is in essence a many-body factor. The many-body nature of bij makes the bond 
energy depend on the local environment of the bond. This feature allows the Brenner potential to predict correct 



IV. CONCLUSIONS 
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geometries and energies for many different carbon structures. The empirical bond-order function used here is written 
as a sum of terms 6,j = [&y _7r + 6^~ 7r ]/2 + bj- where values for the functions ^ _7r and bj~* depend on the local 
coordination and bond angles for atoms i and j , respectively. The first term is a function of the bond angles similar 
to that in the Brenner potential™ while the second term incorporates the third nearest neighbors via a bond-order 
term associated with the dihedral angles and becomes nonzero upon bending of t he grap hene sheets. The values for 
all the parameters used in our calculation for the potentials can be found in Refs. HUH and arc therefore not listed 
here. 

We used the Brenner potential 2 ^ and alternatively the REBO potential by using Large-scale Atomic or Molecular 
Massively Parallel Simulator package 28 ( LAMMPS) for carrying out the molecular dynamic simulations in this work. 
We calculated the pseudo-magnetic fields for the two different potentials along the azimuthal angel at given radius, 
r ~ 40A . Fig. [51 shows the pseudo-magnetic fields B z as a function of along a circular ring of radius r — 4 nm. 
Note that the order of the pseudo-magnetic field is the same in both cases and moreover, we see the same envelop 
behavior at long-wavelength region. However, they exhibit different oscillations modes at short wavelength due to the 
way in which the bond-order ( the dihedral angle) is handled in the REBO potential. 



Appendix B: Simulation methods 



Computer simulations generates very detailed information at the microscopic level and the conversion of this 
information into macroscopic level is the province of statistical mechanics. Molecular dynamics, on the other hand, is 
an important tool to investigate the microscopic behaviors by integrating the motions of particles or particle clusters. 
In the molecular dynamics, the trajectories of atoms are determined by numerically solving Newton's equations of 
motion for a many-body interacting systems, where forces between the particles and their potential energy are defined 
by certain force fields. We used NVT ensemble where the system is isolated from changes in moles (N), volume (V) 
and temperature (T). In NVT, the energy of endothermic and exothermic processes is exchanged with a thermostat. 
A variety of thermostat methods is available to add and remove energy from the boundaries of a MD system in a 
nearly realistic way, approximating the canonical ensemble. On the other hand, a micro-canonical molecular dynamics 
trajectory may be seen as an exchange of potential and kinetic energy, with total energy being conserved. For every 
time step, each particle's position and velocity may be integrated with a method such as Verlet. Given the initial 
positions and velocities, all future positions and velocities can be calculated. If there is a large enough number of 
atoms, statistical temperature can be estimated from the instantaneous temperature, which is found by equating the 
kinetic energy of the system to nksT/2 where n is the number of degrees of freedom of the system. 

We simulated the system at nonzero temperatures by employing a Nose-Hoover thermostat and the time step is 
taken as 0.5 fs. At the beginning of the simulation, we consider a circular graphene ring incorporating atoms which 
are located in the x — y plane. The initial structures are firstly optimized giving the carbon-carbon bond length of 
1.45 A and all samples are initially relaxed at a desired temperature for a duration 1.5 ns. In atomic positions at the 
boundaries, we enforce z = to prevent atomic motion along the z direction and furthermore, after some simulation 
run-times, we set x = y =const. to clamp atoms. 
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FIG. 2: (Color online) Pseudo- magnetic field, given by Eq. [9] in units of tesla induced by shear strains in a circular graphene 
ring. Here, R = 70, Ri = 30Aand 6 = 12°. 
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FIG. 4: (Color online) Projected the gauge filed, for two mentioned cases. Notice that the length value of the vector denotes 
the absolute value of the gauge filed. 




FIG. 5: (Color online) Distribution of the pseudo-magnetic field, B z on graphene ring for two mentioned cases. Green, red and 
black colors correspond to positive, negative and zero values of the pseudo-magnetic field, respectively. For the case one, the 
dark-green(or dark-red) correspond to the strongest pseudo-magnetic filed which is about 1000T and for the case case two, it 
is about 200T 
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FIG. 6: (Color online) Pseudo-magnetic field, B z on a specific direction indicates in the inset figure, left panel, and along 
y-direction passing the origin of the ring, right panel. 
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FIG. 7: (Color online) Pseudo-magnetic field, B z as a function of azimuthal angle for case one at r = 4 nm, left panel and 
r — 6.5nm for case two( right panel). Solid line curve denotes a function which is proportional to cos(3(9). 
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FIG. 8: (Color online) Pseudo-magnetic field, B z as a function of azimuthal angle at T = 50K along a circular ring of radius 4 
nm when the Brenner potential and the REBO potential — are used. Their long- wavelength behavior is very similar to the 
form B z oc cos(3#), however their short wavelength behavior is different due to the way in which the bond-order is handled in 
the REBO potential. The outer boundary of graphene ring, in both cases, is rotated about 9 = 8° . 



